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System and Method for Representing and Processing and 
Modeling Subterranean Surfaces 

FIELD OF THE INVENTION ; 
5 The present invention relates to the field of earth 

models for sxxbterranean surfaces. In particular, the 
invention relates to systems and methods for improved 
representations and processing technic[ues for 
s-ubterranean earth surfaces in earth models used in the 
10 exploration and production of hydrocarbon reservoirs . 

BACKGROUND OF THE INVENTION : 

In the field of processing earth model data for use 
in the extraction of hydrocarbons from the earth, 

15 significant resources have been invested in creating 
functionality and stabilizing the software quality. of 
modeling technology- The efforts have been based on 
faceted representations, and in particular triangulated 
surface methods. However, there are number of key 

20 problems associated with extending triangulated surface 
methods . 

1. Sampling a smooth surface discretely, for 
example with points arranged in a triangle mesh, is 
inherently inefficient- In contrast, smooth surfaces can 

25 be represented as a Taylor series or as an eigen-f -unction 
expansion, e.g., Fourier series of some form. An eigen- 
fiinction expansion can be used to compute an algebraic 
expression to evaluate normal fields, tangent fields, 
etc- This is inherently more compact and efficient than a 

30 2D or 3D sampling lattice with some sort of 
interpolation . 
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2 . The lack of differentiation makes calculation 
of a triangulated surface intersection algorithm 
numerically delicate. 

3. Triangles can be used to' define shape, but 
triangles do not efficiently convey an intuitive sense of 
shape. High resolution and high sampling density make the 
problem more difficult. 

4. Low efficiency triangulation forces application 
developers to either worry about memory management or to 
curb flexibility of data set processing. 

5 . Efficient management of large models means 
localizing change. It is possible in principle to develop 
localization methods using triangulated surface models 
but there are nximerical stability issues seen in 
reference to intersection. 

SUMMARY OF THE INVENTION : 

Thus, it is an object of the present invention to 
provide axi improved system and method for processing data 
used for hydrocarbon extraction. 

Advantageously the invention allows for improved 
memory and CPU efficiency of implicit surface 
construction and editing algorithms. 

According to the invention a method is provided for 
processing data used for hydrocarbon extraction from the 
earth. The method includes the following steps. 
Receiving sampled data representing earth structures . 
Identifying one more symmetry transformation groups from 
the sampled data. Identifying a set of critical points 
from the sampled data. Generating a plurality of 
subdivisions of shapes, the subdivisions together 
representing the earth structures, the generation being 
based at least in part on the set of identified critical 
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points and the symmetiry transformation groups. Processing 
earth model data using the generated subdivision of 
shapes. And, altering activity relating to extraction of 
hydrocarbons from a hydrocarbon reservoir based on the 
5 processed earth model data. 

The identified symmetry transformation group is 
preferably a set of dif f eomorphisms that act on a 
topologically closed and bounded region in space-time 
such that under transformation the region occupies the 

10 same points in space. 

The identified symmetry transformation groups 
preferably correspond to a plurality of shape families, 
each of which includes a set of predicted critical 
points. The subdivisions are preferably generated such 

15 that a shape family is selected from the plurality of 
shape families that corresponds to the identified 
symmetry transformation group. The selection is 
.preferably being based on closeness of correspondence 
between the identified critical points from the sampled 

20 data and the predicted critical points of the selected 
shape family. 

Each shape family preferably has an associated set 
of symmetry transformation group orbits, with each orbit 
being associated with orbit information that specifies 

25 whether the orbit contains a predicted critical point and 
value of the Gaussian curvature of a point in the orbit. 
The orbit information from the set of symmetry 
transformation group prbits associated with the selected 
shape family is preferably applied to the sampled data 

30 thereby generating a unique specification of a shape from 
the selected shape family. Each of the plurality of 
subdivisions of shapes is preferably generated by 
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identifying a part of the uniquely specified shape that 
corresponds to the sampled data. The identified parts 
are assembled, thereby generating a representation of the 
earth structures. 
5 The plurality of subdivisions are preferably 

generated such that the niomber of parameters in each 
subdivision times the number of subdivisions is 
substantially less then would be needed using a faceted 
representation method, and the plurality of subdivisions 

10 is more numerically stable than third order or higher 
representation . 

The invention is also embodied in a system for 
improved extraction of hydrocarbons from the earth, and 
in a computer readable medium capable of causing a 

15 computer system to carry out steps for processing data 
relating to earth structures for the extraction of 
hydrocarbons . 

BRIEF DESCRIPTION OF THE DRAWINGS ; 

20 

Figure 1 shows an open surface and its embedding in 
a closed surface; 

Figure 2 shows a salt mass's steep flanks and 
overhangs ; 

25 Figure 3 shows an example of a 4D representation of 

a field in Turkmanistan; 

Figure 4 is an- image of 'the MacKenzie River Delta; 
Figure 5 shows some combinations of involved 
spherical harmonic polynomials, presented in spherical 
30 polar coordinates; 

Figure 6 illustrates that in most cases a surface 
evolves under mcf to a point; 
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Figures 7a-f forms a series of six images showing 

vpmcf suppressing noises- 
Figure 8 illustrates the orthogonality condition of 

the theorem proposed by Athanassenas; 
5 Figure 9 is an aerial image of part of Big Bend 

National Park, showing the approximation of a plateau to 

a characteristic length scale cone; 

Figure 10 is a side view of a noise cone structure; 
Figure 11 is a satellite image of the Labrador 
10 Trough; 

Figure 12 shows s a sequence of folded sediment on 
the coast of the Gulf of Oman; 

Figure 13 is an image illustrating . progressive 
flattening of an overburden covering a large salt 
15 intrusion; 

Figure 14 is a diagram illustrating the Morse 
theoretical cell deconiposition for a simple configuration 
of a capped and bent cylinder; 

Figures 15-17 show diagrams to aid in the 
20 understanding of the bulls eye construction; 

Figure 18 shows two views of an example of a monkey 
saddle ; 

Figure 19a shows the Reeb graph of a standard torus; 
Figure 19b schematically illustrates a 2D cell 
25 suspension that is induced from the axes and planes of 
symmetry and critical point theory; 

Figure 20 is cross section of the torus shown in 
Figure 19a; 

Figure 21 is a schematic of the shape synopsis 
30 diagram of the torus shown in Figure 19a; 

Figure 22 is a cyclide that is shaded according to 
Gaussian and mean curvature; 
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Figure 23a is a bl-torus with its associated REEB 
diagram; 

Figure 23b is the visual representation of the bi- 
torus shape synopsis diagram; 

Figure 24 is a diagram of the octree with a coarse 
level and leaf level shape index relationship indicated; 

The Figure 25 is a diagram shows part of the French 
model ; 

Figure 26 is an image of the topography of Crater 
Lake, Oregon; 

Figure 27 shows a salt weld in the Gulf of Mexico; 

Figures 28a-c illustrate an example of the misfit 
reduction process; 

Figures 29a-c illustrates an example of blending a 
non-dif f erentiable join of two collars; 

Figure 3 0 shows as a geological exainple a water 
breach as indicated by the white arrow; 

Figure 31 is a NASA Shuttle Mission photograph of 
the Richat Structure in Mauritania; 

Figure 32 shows the natural analog to a conformal 
grid with a proportionally spaced correlation scheme; 

Figure 33 illustrates a non-conf ormal 3D Cartesian 

grid; 

Figure 34 is an image of the Devil's Potholes, South 
Africa; 

Figure 35 is an image of the Yukon River delta; 

Figures 36a and 36b show, for reference, the 
backgroixnd Zechstein Salt and the region in the Zechstein 
where the vsp was acquired; 

Figure 37 show the frame graph that ties the vsp 
region of interest to the Zechstein Salt background; 
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Figures 38a, 38b, 39a and 39b illustrate the 
separation of faulted sediments from unfaulted sediments; 

Figure 40 illustrates a time-lapse seismic 
evolution; 

5 Figure 41 shows the reference structure for spatial 

frames to define a topology graph; 

Figure 42 is a schematic illustration of the 
definition of the variables in the mean curvature 
estimator; 

10 Figure 43 is an image of the regularly sampled input 

surface representing the present example of the top of 
the salt, rendered by drawing a sxxbset of evenly spaced 
inlines and crosslines; 

Figure 44 shows the surface when re- sampled by 
15 interpolating missing sample points; 

Figure 45 shows the re-sampled surface after 25 
iterations of .smoothing; 

Figure 46 shows the shape index map for. the re- 
sampled surface after 25 iterations of smoothing; 
20 Figure 47 shows that there are 33 shapes in the 

example shown in Figure 46; 

Figure 48 is a schematic illustration of a system 
for improved extraction of hydrocarbons from the earth, 
according a preferred embodiment of the invention; 
25 Figure 49 shows further detail of a data processor 

according to preferred embodiments of the invention; 

Figure 50 shows steps in a method for processing 
data used for hydrocarbon extraction from the earth, 
according to preferred embodiments of the invention; and 
30 Figure 51 shows further detail of steps in 

generating an efficient and robust subdivision of shapes, 
according to preferred embodiments of the invention. 
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DETAILED DESCRIPTION OF THE INVENTION : 

While Low-level curvature-based methods as applied 
to implicit surfaces are relatively complicated to 
5 develop, they do not suffer numerical stability problems. 
By comparison, on a smooth Riemannian manifold the 3rd 
derivatives of the square of the signed distance function 
describes the norm of the 2nd Fundamental Form and the 
mean curvature. For triangulated surfaces, this result 
10 is difficult to apply, because numerical evaluation of 
3rd derivatives is not guaranteed to be nximerically 
stable. 

According to the invention, differential geometry 
methods of surface representation will now be described. 

15 Many geoscience phenomena are related to some form of 
fluid flow. If the fluid phenomena under study involve 
surface tension, then mean curvature flow (mcf ) and 
variants such as volume preserving mean curvature flow 
(vpmcf ) are accurate modeling tools'. (For those 

20 unfamiliar with mean curvature flow, imagine the fluid 

front moving normal to each particle of the fluid front.) 
The behavior of mcf and vpmcf are well understood when 
either is applied to a smooth convex surface, star-shaped 
surface, a surface of rotation, or an entire graph. 

25 A reservoir structural framework does not seem to be 

an ideal input to mcf, because noise levels degrade the 
accuracy of analytical approximations and framework 
surfaces in general are not well approximated as convex, 
star-shaped, etc. This perception is erroneous. Our 

30 investigation of curvature-based modeling shows that mcf 
can be used to semi-automate its own noise suppression. 
Given smooth surface data, according to the invention a 
method is provided to decompose that surface into a 
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coimected sum of star- shaped or entire graph or 
axisymmetric surface patches. 

The mathematical fo\indation of mcf is substantial, 
so we seek a unified mathematical description of shape 
5 and its evolution. The concept of a fibre bundle is 

satisfactory. Nakahara presents a formal definition of a 
fibre bundle. Standard examples of fibre bundles are 
vector fields, e.g., velocity fields, and tensor fields, 
e.g., stress fields and elastic fields, evaluated over a 
10 sub-volume. In classical differential geometry, curvature 
properties of surfaces are economically studied in a 
fibre bundle setting. We have found that a 4D fibre 
bundle representation of a reservoir framework is no more 
difficult to write down than is a 3D fibre bxandle 
15 representations. The economy of mathematical 

representation is attractive now from the research view. 
It will be attractive from the engineering view, since 
reuse of concepts limits the amount of technology that 
imist be mastered. 
20 The fibre bundle representation of a surface of 

revolution has the following parts. 

1. A base that is homeomorphic to a loop, 
e.g., a circle, an ellipse, or a polygonal closed 
curve . 

25 2. A typical fibre that is homeomorphic to a 

compact connected part of either a conic section or 
a polygonal line . 

3 . A structure -group that is a group of 
diffeomorphisms that smoothly deforms any instance 

30 of the typical fibre to any other instance. 
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Here are a few shapes that are frequently 
encountered in earth modeling represented as fibre 
bundles . 

a. A torus has a circle base and a circle 

5 typical fibre. Its structure group is the rotation 

group SO ( 3 ) . 

b. A cyclide is a torus, but the structure 
group is the rotation group extended by the scale 
factor diagonal group. Intuitively, it is a skewed 

10 torus. 

c . Any compact planar region R with a 
polygonal boundary 3r has a natural biindle 
structure. To see this, compute the region's 
bounding box B and embed R in B. (Think of B as a 

15 solid.) If there exists a part of 3r that does not 

intersect any of the box faces, then attach a skirt 
S of normal rays to 3. Then S, R, and dB enclose a 
siib- volume V. 

We have found that V is the preferred bundle. The 
20 base of the bundle is 9r and a typical fibre is a 

polygonal line. The region R is formed from the rotation 
of a line segment emanating from the centroid of R and 
joined to the bo\indary 9r. The length of the fibre 
changes instantaneously. The structure group is the 
25 Euclidean group of rigid body transformations. The 
following diagram summarizes this construction. 

Figure 1 shows on the left an open surface, i.e. a 
surface that does not enclose space; on the right, a 
closed region (the large box) is subdivided by the same 
30 surface forming an internal boundary. Figure 2 shows a 
salt mass's steep flanks and overhangs, it also shows an 
example of the cyclide shape in depth imaging. This is an 
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example that challenges existing commercial software. 
WestemGeco commercial processing used to construct this 
velocity model has the common limitation of accepting 
single z-valued height field") data only. According to 
5 the invention we describe below a set of planes in the 
voltame of interest such that a multiple z-value body can 
be subdivided into sections such that each section is 
single valued with respect to one of the planes. In other 
words, each plane parameterizes a section of the 

10 reference multi z-valued surface. Each of these planes 
are equipped with a rotation matrix and translation 
vector so that the application can orient the surface 
section so that it appears to be single valued with 
respect to one of the coordinate planes. (It may be that 

15 the rotated and translated section is normal to the (x,y) 
coordinate plane, which is frequently unacceptable to 
grid-based applications.) 

Material properties of a volume of earth can change 
during an evolutionary process. The time scale of the 

20 evolution can vaary from wall clock to calendar to 
geological record. According to the invention, an 
evolutionary process is represented using a 
generalization of the fibre bundle method that is 
employed for shape. This representation is called an 

25 ""evolutionaDTy process". Here are its components. 

1. The base is a line with period Xb- A 
trajectory on the base serves as a, clock. The 
trajectory's sampling increment has units that are 
30 appropriate to wall clock sampling, calendar 

sampling, or geological record sampling. 
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2. A trajectory has con5)act support. This 
means that the process is defined on a closed at the 
beginning and open at the end boionded interval of 
the helicoid. 

5 3 . Two evolutionary processes Pi and P2 can be 

summed if and only if on a shared inte3rval of the 
base helicoid Pi followed by P2 is identical to P2 
followed by Pi. The sum of Pi and P2 is 
Pl(t) [P2(t) [v] ] , where t is a point in time and v is 

10 a point in the volume domain of interest. 

4. A fibre is a compact path- connected sub- 
volvime. A fibre is associated to each base point. In 
other words the process time stamps every sub- volume 
that it affects. 

15 5. A structure group is a 1-parameter 

subgroup of dif f eomoorphisms that define an 
evolution, e.g., vpmcf. Since the structure group 
consists of dif f eomoarphisms the process must be 
invertible. In particular, the process cannot induce 

20 singularities in material space. It is preferable to 

support singularities, so we allow singularities to 
develop at the end of a process ' s time interval . 

We emphasize that we seek only a geometrical 
25 evolution - not necessarily the true physical evolution - 
that describes part of the structural framework. In the 
following example, we recognize uniformly expemding mean 
curvature .flow creates the shape of the individual 
layers. Each discontinuity in the sediment terminates the 
30 current flow model and is part of the initial conditions 
that define the new flow. The ensemble* of flow problems 
describes the formation, but a more economical 
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representation can be defined if we can treat the entire 
set of restricted flow problems as a single mean 
cuirv^ature flow problem where evolution continues beyond 
the intermediate singularities. We look for an evolution 
5 of the set of initial conditions for the individual flow 
problems, given just the oldest sediment as an initial 
condition. 

We describe a sedimentary sec[uence as a 4D fibre 
10 bundle. 

1. The base of the bundle is a finite length 
piecewise linear curve. 

2. A fibre is an infinitesimal layer of 

15 sediment . 

3. The bxondle's structure group is a set of 
diffeomorphism groups such that each group defines 
an instance of uniformly expanding mean curvature 
flow. Each leg of the bundle base defines a separate 

20 mean curvature flow problem. The transition between 

mcf problems is exactly matched by the discontinuity 
in the sediment. 

Figure 3 shows an example of a 4D representation of 
25 a field in Turkmanistan. Each layer is a distinct mean 
curvature flow where the discontinuity is a curvature 
flow terminator. Working backwards from a flow 
termination, we see that it is much easier to identify 
the flow components in the image. Each flow component is 
30 a uniformly expanding set of solutions to mean curvature 
flow (or volume-preserving mean curvature flow) . Ecker 
shows in his lecture notes that mean curvature flow can 
evolve cracks and holes in a smooth background. See, K. 
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Ecker , Lectures on Regularity for Mean Curvature Flow , 
http : / / web . mathematik . unif reiburg . de /mi /analysis /lehre/WS 
0001/Ecker_WS0001.ps . 

Therefore its usefulness is not limited to modeling 
5 smooth elastic behavior. 

According to the invention, we show how to construct 
a meandering river as a fibre bundle. Figure 4 is an 
image of the MacKenzie River Delta. 

In Figure 4, the river does not intersect itself, 
10 i.e., no oxbow structures are evident. Also, the river 

appears to have constant width. The structure as a fibre 
bundle is clear. 

1. The bundle base is a line segment. 
15 2. The fibre is axxy convenient approximation 

of the river channel in cross section, e.g., a 
trapezoid, 

3 . The structure group is the extended 
Euclidean group, which consists of rigid body motion 
20 plus scaling. 

At each point along the channel we measure the cross 
section and record its position relative to the image 
coordinate axes whose origin is the lower left corner of 
25 the picture. Given the previous position of the cross 
section, we compute the update to the (x,y) plane 
rotation ajad translation. 

We turn now to the question of recognizing the 
elementary shape basis elements in an implicit surface 
•30 representation. Recall that we represent an implicit 
surface as the zero level set of the signed distance 
function (sdf ) . Our constructor solves the signed 
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distance f\inction on a 3D structured grid with tri-linear 
interpolation as a local approximation to sdf , 

By definition tri-linear interpolation T(x,y,z) on a 
grid cell is 

5 

where the coefficients {A^} are defined on the grid 
cell comers . 

10 Each term in T{x,y,z) is an independent 3D spherical 

harmonic polynomial in Cartesian coordinates. For 
convenience we en\imerate these spherical harmonic 
polynomials, using the classical Ymn notation. (See R. 
Baerheim, Coordinate Free Representation of the 

15 Hierarchically Symmetric Tensor of Rank 4 in 

Determination of Symmetiy , Ph.D. thesis. University of 
Utrecht, #159, 1998, Appendix, pg. 141-143). 



Y(m,n) 


Monomial 


Y{0,0) 


1 


Y(0,1) 


z 


Real (¥(1,1) ) 


X 


Iinag(Y(l,l) ) 


y 


Real(Y(2,l) ) 


xz 


Imag(Y(2,l) ) 


yz 


Imag(Y(2,2) ) 


xy 


Iina<i(Y(3,2) ) 


xyz 



20 

Here are the associated symmetry transformations and 
the corresponding isomorphism groups involving these 
spherical harmonic polynomials. 



Term 


SymTnetry generators 


SymTnetry group 


Constant 


Constant 


so (3) • 


X 


[X] - [x] 


S0(2) 


Y 


[y] ~ [y] 


SO (2) 
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z 


tz] - [z] 


SO (2) 


X + Y 


[X, y] ~ [y, X] 


Z/2Z 


X + Y + XY 


Y + Z 


[y, z] ~ [z, y] 


Z/2Z 


Y + Z + YZ 


X + Z 


[x, z] ~ [z, x] 


Z/2Z 


X + Z + XZ 


X + Y + Z 


[x," y, z] ~ [y, 

2, Xj 

[x, y] lYr x] 


Tetrahedron 
group of order 
12 


XYZ 


[x, y, z] - [y, 

Z / XJ 

[y/ z] - [-y, -z] 
[X, y] [y, x] 


Symmetric (4) 


XY 


Lx, yj — L-x, — yj 


Z/2Z e Z/2Z 


[x, y] ~ [y, x] 


YZ 


[y, z] ~ [-y, -z] 


[y, z] - [z, yJ 


XZ 


[x, z] - [-X -z] 


[x, z] - [z, x] 


XY + XZ 


[x, y, z] [-X, 

-y, -z] 
[y, z] ~ [z, y] 


Z/2Z e Z/2Z 


XZ + YZ 


[x, y, z] - [-X, 

-y, -z] 
[x, y] - [y, x] 


XY + YZ 


[x, y, z] - [-X, 

-y/ -z] 
[x, z] [z, x] 


XY + YZ + ZX 


[x, y, z] - [-X, 

-y. -z] 
[x, y, z] - [y, 

z, x] 
[x, y] [y, x] 


Symmetric (4) 



Symmet3:Y analysis is a form of spectral analysis 
applied to the discrete spectrum that is associated with 
spherical harmonic polynomial expansions. Our estimates 
5 of tri-linear interpolation coefficients are noisy, so we 
need a threshold for dismissing spectral lines. 

Figure 5 shows some combinations of involved 
spherical harmonic polynomials, presented in spherical 
polar coordinates . 
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We are interested in conic section fibre bundle 
shapes. Here is the correspondence of shape to symmetry 
group . 



Symmetry group 


Candidate shapes 


SO(3) 


Sphere 


0(2) 


Cylinder, spherical torus, 
elliptical torus, hvperboloid 


S0(2) 


Cone, paraboloid 


Symmetric (4) 


Ciibe 


Dihedral (4) 


Ellipsoid 


Dihedral (3) 


Triancrular prism 




Rectangular prism 


Z. 


Cvclide 




Noisv data 



According to the invention, mean curvature flow and 
volume-preserving mean curvature flow will now be 
discussed in further detail. We define a few terms that 
appear frequently herein. Given a 2D manifold M and a 
10 point p E M, let {pi, Pa} be a local coordinate system for 
a region containing p and let n be an outward pointing 
normal at p. Finally the Euclidean inner product of 
vector Ui and Vj is denoted by <Ui, Vj>. 

Definitions 

/dM dM\ 

The First Fundamental Form (1 FF) is = V^p'''dP'/' 
•^^ The inverse of the 1" FF is denoted by 

T'he Second Fundamental Form ^2"" FF) is A = [ft^] = -( gj^.^y .»/• 
The Weingarten map is W/ = g^^h^,. (Einstein notation used.) 

The eigenvalues of the Weingarten map are the 
principal curvatures. The trace of the map is the meMi 
curvature , the determinant of the map is the Gaussian 
20 curvature . The norm of the 2^ FF jAl^ is defined as the 
sum of the squares of the principal curvatures. 
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We define mean curvature flow and volume-preserving 
mean curvature flow. 

Notation 

(M^) is a family of evolving smooth 2D manifolds such that 
M = is given. 

N(x, t) is the normal at x e M^. 

H(x,t) is the mean curvature evaluated at x 8 M,.. 
h(t) is the average value of H(x,t) on M^. 

Mean curvature flow (mcf ) is defined as the solution 
to the initial and boundary value problem 

^ = -H(xj)'Nixj). [4] 

VolTime-preserving mean curvature flow (vpmcf ) is 
defined as the solution to the initial and boxindaxY value 

problem —L = -(^H(xa)-h(t))'Nix,t). [5] 
ot 

Figure 6 illustrates that, except in ideal 
circiimstances (when the input is a smooth closed convex 
region or a graph) , a surface evolves \mder mcf to a 
point. Technically, mcf develops a singularity in finite 
time. See, K. Museth, D. Breen, R. Whitaker, A. Barr, 
Level Set Surface Editing Operators , SIGGRAPH 2002. 

A straightforward way to prevent annihilation of 
this shape is to stop the mcf after some number of time 
steps, inspect the results, and maybe resume the process. 
This is not convenient in a production environment. In 
fact, the manner in which uncontrolled mcf annihilates a 
shape enables us to attach a recognition procedure that 
mcf can call to decide when to ask the human for 
permission to resume the figure's evolution. 
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Figures 7a-f forms a series of six images showing 
vpmcf suppressing noise. See, A. Kuprat, A. Khamayseh, D. 
George, L. Larkey, Volume Conserving for Piecewi se Linear 
Curves, Surfaces, and Triple Lines / Journal of 
5 Computational Physics, 172 (2001), pg. 98-118. In Figure 
7a, the southern hemisphere is corrupted with a 
significant amount of Gaussian noise. By Figure 7c, it 
makes sense to consider how much additional noise, if 
any, must be removed before the smoothed result is an 

10 acceptable approximation . 

We are interested in surfaces of revolution. The 
following result has been obtained regarding the 
behaviour of surfaces of rotation \inder vpmcf. See, M. 
Athanassenas, Volume-preserving mean curvat ure flow of 

15 rotationally symmetric surfaces , Comment. Math. Helv. 
72(1997), pg. 52-66. 

Theorem (Athanassenas) 
2 0 Let Mo be a smooth rotationally symmetric surface 

enclosing a sub-vol\ime V. Let S be a slab of thickness X 
> 0 that is bounded by the z = 0 and z = x height field 
planes. Suppose that Mo satisfies the following two 
conditions . 

25 

a. Mo has a line of intersection in the z = 0 and z = 
T height field planes. 

b. The end points of each line of intersection 
between M^ and S meet an end of M^ at a right 

30 angle. 

If the voliame of V is greater than x times the area 
of M, then Mt evolves \inder vpmcf to a cylinder. 
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As an example of the condition on the lower bound on 
volvime, consider the case of a right cylinder of height x 
and radius r. This cylinder has volxame V equal to TCr^x, 
while its surface area M is equal to 2nxx. Hence the lower 
5 boxuid on the volume is satisfied exactly when Tir^x > 

27tri;^, i.e, r > 2t. Heuristically, the bound is satisfied 
for ^"sccuat" cylinders. We observe that a sphere cannot 
satisfy the volvme to area relationship. To do so would 

imply that there exists r > 0 such that — ^r^/4^r^ =-^>2r . 

3 3 

10 Figure 8 illustrates the orthogonality condition of 

the theorem proposed by Athanassenas . In this diagram we 
show two vertical cross-sections. In the figure on the 
left the intersection of the figure with the top and 
bottom planes must satisfy the right angle hypothesis. 

15 The curvature flow takes care of the irregular vertical 
surfaces either end of the figure. In finite time vpmcf 
generates from the figure on the left the figure on the 
right . 

Volume preservation is essential for earth model 
20 applications, so we prefer vpmcf to ordinary mcf . Later 
in the description, we will use this theorem to reduce 
the discrepancy between an idealized representation of 
shape and a sampled data surface. 

Another important class of smooth surfaces are those 
25 that are star- shaped. A surface is star- shaped if there 
exists a point P on the surface such that a line segment 
PQ that is entirely contained in the surface can join 
every other point Q in the surface. 

It has been shown that star-shaped closed smooth 
30 surfaces are stable under mcf. See, K. Smocyzk, Star- 
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20 



25 



30 



shaped hvpersur faces and the mean curvature flow, 
Manuscripta Math. 95 (1998), pg. 225-236. 

We show now that vptncf suppresses additive high 
frequency harmonics before decaying the imderlying low 
frequency shape signal. Vpmcf enjoys the property that 
surface area is always decreasing in time. 



10 We frequently need to compute mean and Gaussian 

curvature for a single valued surface over a plane, i.e., 
a graph. Suppose that S = S(x,y). Then the formulae for 
mean curvature H and Gaussian curvature K are as follows 



2 

nr — "^">y and 

15 ^ - 77"; 2 . . 2n2 



We model a surface as a collection of patches, where 
each patch is taken from a surface of revolution, say S = 
(${v)*cosu, <|>(v)*sinu, V(v)), such that v e (a,b) , u 
e(0,27t), and \|r(v) ?t 0. The formulae for mean curvature H 
and Gaussian curvature K for a surface of revolution is 
as follows. 



2 



For the proofs, see, M. do Carmo, Differential 
Geometry of Curves and Surfaces, Prentice Hall, 1976, 

pages 162-163 . 

We need to understand how mcf behaves when applied 
to a noisy surface S, where S- is either single valued 
over a plane or is a surface of revolution. The above 
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formulae imply that if S is an entire graph, then the 
harmonic number N scales the mean curvature. Similarly, 
if S is a surface of revolution, then scales the mean 
curvature. (Spherical harmonic polynomials have a cosine 
5 factor, so these estimates remain valid under the natural 
heat equation eigen- function expansion.) Consequently, 
the mean curvature integral involves a scale factor of 
either or N^. Therefore the rate of change of surface 
area becomes very negative as N increases . We conclude 

10 that vpmcf eliminates high frequency harmonics during the 
initial stages of the evolution ahead of low frequency 
harmonics that form the ideal shape. A maximum or miniimm 
in the ideal shape will. appear to move as vpmcf 
eliminates corrupting noise, but the critical point's 

15 location will stabilize as the noise disappears. 

Importantly, when a critical point's location stabilizes, 
it does not move later in the evolution unless vpmcf 
begins to decay the shape. We remark that mcf can be 
turned off locally by resetting the mean curvature to 

20 zero. 

We turn now to the question of how to decide when 
noise suppression turns into shape decay. Refer again to 
Figures 7a-f , which show the smoothing of a noisy sphere 
\inder vpmcf. Athanassenas' s theorem does not apply, since 
25 the sphere fails the vol\ime bound assumption. However, it 
still makes sense to apply Kuprat et al.'s vpmcf 
procedure to the shape. When we do this we get Figures 
7c-f . 

It is reasonable to say that noise suppression 
30 reduces Figures 7a-b to Figures 7c-d, It is harder to say 
if image noise suppression or shape decay accovmts for 
the transformation to Figures 7e-f . Therefore, we want 
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vpmcf to proceed without manual intervention to eliminate 
noise but to seek guidance when the flow causes the 
\inder lying shape to decay. Here is a way to monitor the 
elimination of noise and detect decay of the critical 
5 points in an underlying shape that relies on vpmcf. 

1. Given a characteristic length, we embed a series of 
conical shapes in the surface. The base of the cone 
is either the outermost ID boundary of the surface 

10 or it is an inflection point curve surrounding the 

patch that contains the maximum or minimum. 
2 . We recognize minima and maxima in its 2D enclosure 
using a standard nested bounding box algorithm. If 
the surface is implicit, then we can localize the 

15 region by testing the implicit surface's 

triangulated parent. 

3. Use this 2D enclosure as the top of the (noisy) 
cone. Note that we cannot assume that each maximum 
or miniimim is topologically isolated, e.g., a ridge 

20 might exist. 

4. Apply any mean curvature flow procedure, e.g., D3b, 
to the surface and monitor the effect on the cone's 
slant height. Mean curvature flow theory says that 
the slant height decreases as the noise is removed 

25 (it becomes straighter) . However, when the mean 

curvature flow begins to destroy the intrinsic shape 
of the surface in the vicinity of the cone, then the 
slant height will decrease beyond a threshold. 
Therefore, the test for adequate smoothing is to 

30 compare the ratio of slant height before smoothing 

to the result obtained after each iteration. 

Figure 9 is an aerial image of part of Big Bend 
National Park, showing the approximation of a plateau to 
35 a characteristic length scale cone. It is difficult to 
precisely locate the maxima on the plateau, but it is 
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easy to enclose the region containing the maxima in a 
tight loop- It is unimportant that the cone does not have 
a circular cross section. 

Figure 10 is a side view of a noise cone structure. 
5 When monitoring noise suppression on the surface of a 
geobody, e.g., it is useful to have a noise-monitoring 
device for thin undulating cross section. We replace a 
torus by a cone. We use the cone to monitor shear 
stretching and erosion of an interface. A clear instance 
10 of this phenomenon is Figure 7c. 

According to the invention, mean curvature flow and 
framework I/O will now be discussed in further detail. 
Singularities mark an end to the smooth evolution of a 

15 shape under mcf . A singularity is frequently easier to 
identify on an image than is the interface that 
' corresponds to the precise start of the flow. The flow 
imposes a natural partition of the framework. Each region 
in the partition is a self-contained expression of mcf. 

20 When we think about sending and receiving an update to a 
framework, we prefer to send and receive a mcf problem 
with boiindary conditions rather than an opacjue byte 
stream. We specify the solution form - whether the flow 
unifo2™iy expands or contracts the solution or smoothes 

25 the initial data -and provide beginning and end surface 
data. We describe the intermediate surfaces by recording 
the time states that correspond to the intermediate 
surfaces. We also send critical point data for each 
intermediate surface. Interested applications rim the 

30 vpmcf script and reconstruct the update locally. We trade 
fast CPU for less fast I/O. 
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Here are three illustrations of this idea. Figure 11 
is a NASA satellite image of the Labrador Trough. In this 
image we notice that the sediment resembles a 
longitudinal cross section of a brain with an attached 
5 spinal cord. Singularities separate "brain tissue" from 
the stem of the « spinal cord". A singularity also 
s\ibdivides the "cerebellum" . 

Figure 12 shows s a sequence of folded sediment on . 
the coast of the Gulf of Oman. We model this as a 
10 uniformly expanding (seen right to left) solution to mcf . 
Ecker proved that a uniformly expanding solution to mcf 
is given by the equation M , = ^/T • M , , provided that the 
initial fold is approximately a cone. 

Figure 13 is an image illustrating progressive 
15 flattening of an overburden covering a large salt 

intrusion. We notice that the fault block that occupies 
the left half of the image (indicated by the lower white 
arrow) displays a gradual flattening in space and time as 
the top of the salt's intrusion is progressively reduced 
20 as the sediment was deposited on top of the salt 
(indicated by the upper white arrow) . 

According to the invention, a preferred technique of 
2D parameterization will now be discussed in further 

25 detail. Many surface editing operations are more 

efficient when the operator can access the surface's 
parameterisation. Height field data uses the (x,y) co- 
ordinate plane as the parameter space and coordinate 
projection as the parameter space mapping. We explain how 

30 to obtain a parameterisation of a smooth surface S with 
no assuits>tions regarding the orientation of the surface 
relative to the (x,y) co-ordinate plane. 
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We claim that there exists an invertible projection 
of S onto a cylinder that has 0, 1, or 2 end caps - 

1. Given any two closed curves and on S. 
We partially order the two cxirves^ say < C^, if 
5 there exists a 2D region on S such that dn^ = 

and CiCR, We call the 2D region of S that is 
contained between \ the collar that is bounded 
by Cj and 

2 • From the classification of smooth surfaces 
10 in we know that S has Euler Characteristic X(S) 

eciual to 0, 1, or 2 . The Euler Poincare Theorem says 
that %(S) = # (maxima) + # (minima) - # (saddle 
points) - Note that we use height as our Morse 
f\inction, so the term critical point" is a 
15 contraction for ^'height field critical point". 

3 . We will prove our claim by induction on 
the number of saddle points on S. 

4 . Suppose that the number of saddle points 
on S is zero. 

20 a. If x(S) = 2, then Reeb's Theorem says 

that S is dif feomorphic to a round sphere, 
which is dif feomorphic to a cylinder with two 
end caps . 

b. If %(S) = 0, then S must be 

25 dif feomorphic to a plane which is dif feomorphic 

to a cylinder with one end cap. 

c. If x(S) = 1 then S has either 1 
maximum or 1 minimum but not both. 

d. Construct a 1*^ order Taylor series 
30 with center equal to the critical point and 

quadratic order error term. 

e. Notice that the gradient term is 
trivial, since we expand about a minimum or 
maximum. 
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f . We apply a translation to the surface 
so that expansion point is the origin and has 
height equal to 0 . 

g. We conclude that after a suitable 

5 translation and (x,y) plane rotation a 2 term 

Taylor expansion about the critical point with 
a quadratic error estimator siitvplifies to a 
pure quadratic t-erm. 

h. Given 8 > 0, we fit the Taylor 

10 expansion such that the error estimate is less 

than e. 

i. Then the Taylor series esqpansion 
within Co is a cap/cup. We know that a cap/cup 
is diffeomorphic to a plane and that a plane is 

15 diffeomorphic to cylinder with one end cap. 

Furthermore, the open end of the cylinder has 
boundary C^. 

5. Suppose that the number of saddle points 
20 on S equals 1. 

a. If X(S) = 2, then S possesses 
, either, (i) 2 maxima and 1 minimum or (ii) 2 
minima and 1 maximum or (iii) 3 maxima or (iv) 
3 minima. 

25 b. Case (iii) is impossible, because 

Morse Theory says that a maximtom contributes a 
2D cell to the shape. If case (iii) were true, 
then 

c. Case (iv) is impossible, because 
30 Morse Theory says that a minimijm contributes a 

OD cell to the shape. If case (iv) were true, 
then S contains no 2D cells. 

Figure 14 is a diagram illustrating the Morse 
35 theoretical cell decomposition for a simple configuration 
of a capped and bent cylinder. 
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Morse Theory explains the cell decomposition of the 
shape. (See M. Hirsch, Differential Topology^. Springer 
Verlag, 1976, pg. 156-164.) 

5 d. This shape has a saddle point at 

location C. 

i , Morse Theory says that the 

surface contains a ID cell that is 
attached to location C. 
10 ii. We observe that this ID cell is 

not a loop. Specifically, the point at 
location D is not part of the ID cell. 



15 D. 



e. This shape has a miniimim at location 

i. Morse Theory says that the shape 

contains a OD cell that is attached to 
location D. 

20 f • This shape contains two maxima at 

locations A and B. 

i. Morse Theory says that the shape 

contains two 2D cells, and that each 2D 
cell is attached to a maximum. 

25 ii. The boundary between the two 2D 

cells is the union of the OD cell that is 
attached to location D and the ID cell 
that is attached to location C. 
iii . We decompose each 2D cell into 

30 the union of a Taylor series cap expainded 

about the maximum and a collar that is 
bo\mde4 by the loop suspended between 
locations C and D and the Taylor series 
circle of convergence. 



35 



g. Suppose that X(S) = 1. The Euler- 
Poincare formula says that the surface contains 
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either (i) 1 maximum and 1 minimum or (ii) 2 
maxima or (iii) 2 minima. 

i. An exait^le configuration for 

case (i) is the upper half of a torus . 
5 ii. An example configuration for 

case (ii) is a flat surface with two 
hills. 

iii. An example configuration for 

case (iii) is a flat surface with two 
10 valleys . 

h. Finally suppose that X(S) = 0. Then S 
contains either 1 minimiam or 1 maximum, but not 
both. An example of this configuration is a 
15 single mountain or sinkhole attached to a 

plane. 

6. We have specified the parametric mapping 
when S contains either 0 or 1 saddle point. 



20 



25 



30 



35 



7 . Now suppose that S contains N-1 saddle 
points. Choose any saddle point and pass a height 
field plane through it, subdividing S in two. At 
least one of the two parts contains fewer saddle 
points than does S. Define S* to be this part of S. 

a. By induction we know that S* is 
diffeomorphic to either a cylinder with a hole 
in its lateral surface, or a cylinder that has 
0 or 1 end caps . 

b. Closing the hole in S* corresponds to 
either closing a hole in the walls of a 
cylinder or to attaching an end cap to the 
cylindrical surface that we created in xyz 
space. 

Figures 15-17 show diagrams to aid in the 
understanding of the bulls eye construction. In Figure 
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16, we invertibly map the possibly patched A to the 
planar unfolding of its interior facetted cylinder. 
Figures 17a and 17b illustrates how to perform the bulls 
eye construction in the region of a saddle point. The 
5 saddle point is indicated with the dot 10 in Figures 17a 
and 17b This process is preferably repeated for every 
saddle point of the shape. 

According to the invention, tri- linear 

10 interpolation, implicit surfaces, and critical points 

will now be described in further detail. In particular 
we herein discuss the analysis of the tri-linear 
interpolation approximation to the signed distance 
function on a 3D structured grid. We derive necessary two 

15 conditions that the tri-linear interpolator must satisfy 
if the grid cell contains a height field critical point. 
One condition specifies a relationship between the tri- 
linear coefficients and the critical value. The other 
condition establishes a second relationship among four of 

20 the tri-linear coefficients. Taken together, we show that 
a height field critical point for an implicit surface 
must be embedded in the grid cell faces and can never be 
found in the interior of the grid cell unless the 
implicit surface is a plane. 

25 Let G denote a 3D structured grid with typical grid 

cell C . We assume that an implicit surface S is contained 
in G and in particular that C ns 0. We denote the tri- 
linear interpolation function on C by a and we assume 
that C is small enough that a is a good approximation to 

30 the signed distance function that implicitly defined S. 
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We begin with the definition of tri-linear 
interpolation 

crix,y,z) = Aix,y)+ Bix,y)z d) 
where 

B (X, y) = bo + b^x + b^y + b^xy. 

For later reference, we note that A and B are 
harmonic , i.e., 

VU = V'J5 = 0 

We are interested in the zeroth level of a. If o = 
0, then either B = 0 or B 9t 0 . We suppose that B 0. 
Then 

, A(j^-y) (4) 

Bix,yy 

We know that height field critical points are found 

when 



15 



25 



30 



dx dy 

This condition is satisfied when 



B 4^.- A ^= 0 

d X ox 

d y d y 

Rewriting (6) 



(6) 
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a A ^ A 



A 




d X 




d 


y 


B 




d B 




d 


B 



(7) 



d X d y 



15 



Using (7) 



3 aM_) = 



d y d X d X 

(_ ) + (B a3 - A 63 ) = (8) 

d y ox ax o y 

(S - A &3 ) = 0 . 

The Hessian matrix of 2"^"^ partial derivatives 
10 evaluated at the critical point is 

0 z _ dx dx ax _ A 

dx' B' 

o z _ dy dy dy ^ q 

dy' B' 



(9) 



d^Z dxdy dxdy Ba^ - Ab^ ^ 

dxdy B' B^ 



We conclude that the Hessian matrix of 2"^^ partial 
derivatives is identically zero. We are led to consider 
what kinds of shape have the property that in a local 
20 neighborhood of a critical point the matrix of second 
partial derivatives is identically zero. 

We have enough information to compute the Gaussian 
curvature of z = z(x,y) . Recall that for a graph z = 
z(X/y) that its Gaussian curvature K(x,y) is 
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20 



25 



^(^.3,) = J£13ZL_3||>_ (10) 
9;c oy 



We notice that the signum of K(x,y) depends only on 
5 its numerator. Siibstituting (3) into the numerator of 
(9), we discover that the numerator of K(x,y) is 

10 Hence the Gaussian curvature K(x,y) < 0 everywhere. 

As an illustration we consider the "monkey saddle" 
M(x,y), which coincidentally is also a spherical harmonic 
polynomial M(x,y) = Re(Y33) = -x^ " 3xy^ The monkey saddle 
has a single critical point, which is the origin (0,0). 

15 The matrix D of 2^'^ partial derivatives is 



D = 



6x -6y 
-6y -6x 



(12) 



Sxabstituting (x,y) = (0,0), we find that the 2^ 
partial derivatives matrix is identically 0. We conclude 
that the origin is a degenerate saddle point for M(x,y) . 

Figure 18 shows two views of an example of a monkey 
saddle. In particular the monkey saddle is M(x,y) given 
by the equation above . The image on- the left is shaded 
according to mean curvature, while the image on the right 
is shaded according to Gaussian curvature. We have placed 
a white dotted circle around the critical point. 

The critical point is a degenerate saddle point. We 
see this by substituting the samples defined by the white 
30 circle and observing the pattern of +/- signs. As an 

exan5)le, let 0 < a < b. Then ((+/-) a, { + /-)b) produces the 
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5 



pattern (-,-)/ (-/+)/ (-/-)# and (-,+) clockwise from 
Quadrant I . 

The Gaussian curvature of the monkey saddle is 
(Gray, pages 382-383) 



(See, A. Gray, Modem Differential Geometjcy of 
Curves and Surfaces with Mathematica 2nd Edition^ CRC 
10 Press, 1998, pages 382-383) . 

Summarizing the analysis for B 5^ 0, we have shown 
the following. 

1. A grid cell contains at most one height field 
15 critical point. 

2. If a grid cell does contain a critical point, then 
the critical point's coordinates are uniquely 
determined from the tri-linear interpolator's 
coefficients. We use this as a quick test for 

20 assessing the existence of a critical point in a 

grid cell. 

3 . The matrix of 2**^ partial derivatives is 

identically zero when evaluated at the critical 
point . 

25 4. The Gaussian curvature K(x,y) < 0 and equals zero 

exactly at the critical point. 
5. z = z(x,y) is harmonic, so by the Maximxam 

Principle for harmonic functions on compact sets, 

we know that the minimxom and maximum value for z 
30 on a grid cell will be found on the grid cell 

comers . 
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It has been remarked that the behavior of the tri- 
linear interpolation function is determined in the region 

dz Bz j^z _ 

5 of a critical point by the behavior of ' 

G. Weber, G. Scheuermann, H. Hagen, B. Hamann, Exploring 
Scalar Fields Using Critical Isovalues, http: 
/ /graphics . cs .ucdavis . edu/-hamann/WeberScheuermannHagenHa 
mann2002.pdf ("Weber"). We agree, which is Morse Theory 

10 should be applied. 

Weber shows that a tri- linear interpolation f\anction 
can have a critical point on the edge of a grid cell if 
and only if the tri-linear interpolation function is 
constant along the edge. Weber does not assume that the 

15 tri-linear interpolation approximates a signed distance 
function. Weber also obtains a simple test for the 
existence of a maximum at grid cell vertex (0,0,0) by 
looking at the value of the function at the tetrahedral 
comers (1,0,0), (0,1,0), and (0,0,1). 

20 In our situation Weber's condition says that a© > max 

(ai, as, bo) irt5)lies that (0,0,0) is a maximum of T(x,y,.z) 
= A(x,y) + B(x,y)*z. We note that T(0,0,0) = ao and that 



11 
dx 

dT_ 

By 

dz 



= 0 

(0,0.0) 

= 0 (14) 

(0.0.0) 

= J5(0,0) = bo 

(0,0.0) 



25 



Checking 2"* mixed partials 
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10 



20, 



d^T d-A d^B - 

J— r — Q 



3jc93; d X d y d x d y 



« 3 



a X a z ax 
a d B 



d y d z d y 



^ 1 + ^ 3 y (15) 



Therefore the first order Taylor series expansion of 
the tri-linear interpolation function about (0,0,0) 
inside a grid cell cube is given by the expression 

T (x, y , z) = ^0 + ^oZ + + i^3y)^z + (br, + b^x)yz 

= ^0 + b^z + b^xz + b^yz + b^xyz + b^xyz (i6) 
= ^0 + Bix, y)z + b^xyz 



This says that in a local neighborhood of (0,0,0) 
that the tri-linear interpolation function is a harmonic 
polynomial and therefore the Maximum Principle applies. 
15 We conclude that on the tetrahedron formed by the four 
grid cell vertices that maximum occurs at the grid cell 
vertex ao or bo. We invoke Weber's assumption and conclude 
that the maximum occurs at ao- It is not necessary to 
assume the other inequalities in Weber's assiimption. 



We consider now the case that a = 0 and B = 0 . We 
begin by obseirving that a = 0 implies that B = 0 if and 
only if A = 0. Therefore, 
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an + a.x bQ + b^x 

y = ^ - = - .V , • (17) 

a2 + (^x b2+o^x 

Again we look for solutions to the zero gradient 
equation . 

5 

dy_ ^ _ a^a^ - a^a^ ^ q implies. that a,a^ - a^a^ = 0 
dx (aj + a^x) 

and (18) 

dy_^_b,K^-b^b^^ =0 implies that b^b^-b^h = 0 
dx (&2 + ^z^y 

Again we obtain a relationship among the tri-linear 
interpolator coefficients that must be satisfied in order 
10 that a grid cell face contain a critical point on a curve 
running through the crrid cell face. 

Suppose that a grid cell face contains a critical . 

point of a curve. Again the second derivative _ o, so 

15 we characterize the nature of the critical point by 
sampling the gradient on a tight neighborhood of the 
critical point. 

We conclude that there exists at most 1 critical 
20 point on a curve running through a grid cell face and 
that the (x, y) coordinates of a critical point are 
uniquely determined in terms of the trir linear 
interpolator ' s coefficients . 

25 According to the invention shape index and shape 

identification techniques will now be described in 
further detail. According to a preferred embodiment, 
shape identification depends on a dimension - independent 
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measure of the principal curvatiires at a point known as 
the Koenderink and van Doom shape index si . 



10 



2 ^ + y 
SI = — arctan — 

where K y > /r,. 



— = — arctan — j= 

2 ^ -Jh 



H 



K 



Shovm below is a chart of the shape index map's 
range [-1.0, +1.0], which goes from most concave to most 
convex. We note that a shape index value of zero 
corresponds to a zero mean curvature surface, which for 
the case of a compact surface, equates to a catenoid or a 
compact planar region. 



Local shape 


si interval 


spherical cup 


[-1.-7/8) 


trough 


[-7/8,-5/8) 


rut 


[-5/8,-3/8) 


saddle rut 


[-3/8,-1/8) 


saddle 


[-1/8,1/8) 



Local shape 


si interval 


saddle ridge 


[1/8,3/8) 


ridge 


[3/8,5/8) 


dome 


[5/8,7/8) 


spherical cap 


[7/8,1] 



Cantzler et al . have computed a correspondence 
between shape attributes defined by a Gaussian and mean 
15 curvature value pair and the shape index map range. See, 
H. Cantzler, R. Fisher, Comparison of HK and SC curvature 
description methods , 

http : / / www. dai . ed . ac . uk/homes/rbf /hcSdim.ps . gz . See, the 
following tables. 
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K<0 


K=0 


K>0 


H<0 


Saddle Wley 


Concave. 


Concave 




(SvHy) 


Cylinder (-C)') 


Ellipsoid (-Q) 


H=0 


Miniinal 


Plane 


Imposable 




(MHy) 


(PI) 




H>0 


Saddle Ridge 


Convex 


Cbm-ex 




(SrHy) 


Cylinder (+Cy) 


Ellipsoid (-(-El) 



Shape 


Index range 


Concave Ell%>soid (-El) 


S €[-1.-5/8) 


Concave Cyiinder<-Cy) 


S € [-5/8,-3,'8) 


Hyperboioid(Hy) 


S € [-3/83/8) 


Com'sx Cylinda- (+Cy) 


S€[3/S^./8) 


Convex Ellipsoid C+EI) 


S € [5/843 



As an example, we compute the shape index of a 



torus 



10 



15 



20 



The toirus's parameterization is 

((a + 6 • cos v) cos m, (a + b- cos v) sin m, & • sin w), 
where we assume that a > b > 0 . 

cosv 



and =■"•7- 
a+fccosv b 



25 



The principal curvatures are 
Therefore — ^ = -(1+ — cosv). 

This q[uotient covers the range [-3.0, +1,0]; hence 
the shape index covers the interval [-0.80, +0.50], 
which corresponds to '^trough'' to ^^ridge". 

Certain principal curvature pairs are distinguished. 

1. If K2 = iq, then the shape index is +/- 1.0, as 
the signiamCKj) = -/+1. 

2. If 0 = > iq, then the shape index is 

3. If Kj > iq= 0, then the shape index is Vi. 

4. If = -iq, then the shape index is 0. 

5. The shape index of a maximum is 1-0. 

6. The shape index of a minimiam is -1.0. 

7. If the shape index is discretized in increments 
of size As , then the principal curvatures will 
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be discretized according to the following 
f orntulae 

Shape index expresses a relationship between 
principal curvatures. For example, we expect the shape 
index in the region of an inflection point is 
approximately 1/8. Plugging into the shape index 
definition, we discover that 



10 = /Tj 



tan(— 

2 8 

tan( — ) + 1 



= -0.99023 • K, 



The shape index enables a h\iman to express a 
threshold change in curvature in a- dimension- independent 
manner. This correspondence is fundamental to the 
15 robustness of our shape identification scheme. 

Measurements are always tainted with noise. Therefore it 
preferable to identify intervals rather than point values 
for attribute correspondence. 

Definitions 

Let A,B eS and let s^^Sg be the shape index at A,B. 

20 By definition tan(— • ^y"^) = and t2Ln{^' s^) = 0^ .where 

2 2 

A 



qa _ ^1 + ^ y^here kJ" and kJ" are principal curvatures at A. 
Likewise for 6^ . 
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This formula relates the shape index principal 
curvature quotient in shape index increments. For example 

if is' -^^) =i r/ten tan(^) = 0.0034 = ^f^^^ ■ 

Holding 0"" fixed, then t&ni^)- (1-0^0*) = 0' -0^implies that 

16 

S , tan(— ) + g ^ 0.0034 + ^'' 
16 

Now we can utilize the formula relating the value of 
a shape index to the ratio of the constituent principal 
curvature . 

10 

According to a preferred embodiment of the 
invention, techniques for curvature-adaptive sampling of 
a smooth surface will now be described in further detail - 

15 Due to popularity of 2D FFT methods, it is common 

practice to sample a surface based on a fixed spatial 
step size. For smooth data, we know that a low order 
Taylor polynomial expansion can supply sample data to 
whatever a priori precision and at whatever spacing is 

20 acceptable. There is no need to store a large array when 
an algebraic expression can be evaluated on demand. Our 
method for curvature-based sampling is straightforward. 

1. We are given a cloud of points. Using Hoppe's 
25 tangent plane fitting procedure (Hoppe, pages 

22-25) , we loop through the input point cloud 
fitting tangent planes. 
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2. At each tangent plane, we determine if the 
parent point is a height field critical point. 
We record the set of all critical points {C,,} . 
We use the symbol T(Cj^) to denote the tangent 

5 plane at C,^ 

3 . At each critical point C^, we define a 
neighborhood N{C^) about in the point cloud 
such that the tangent plane T(C,,) is a 
parameterization of the triangulation of the 

10 points in N(C,^) . Abusing notation slightly, we 

use N{C,^) to denote the triangulation of the 
point cloud. 

4. Using equations (4.7^) and (4.72) estimate 
both principal curvatures cind the associated 

15 shape index on N(C^) . 

5 . Given a shape index increment , we apply the 
analysis developed in Section 8 to partition 
N(C^) into shape index equivalence classes. We 
chose to seed the shape inde?c equivalence 

20 classes with critical points so that there is 

the least ambiguity in the shape index 
assignment. 

6. If the critical point neighborhoods do not 
cover completely the implicit surface, then 

25 search the remaining surface. for regular points 

that have an iinambiguous shape index interval 
assignment. 

7 . Bin the remaining unassigned sample points to 
new shape index equivalence classes. 

30 8. Assign each shape index equivalence class to a 

distinct grid cell. 
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Done, 



This algorithm uses just curvature information, 
since we apply it to the unorganized input point cloud. 
We can improve the N(Ck) partition if we know a priori 
that the shape identification procedure described above 
has been applied. Using the information generated during 
shape identification, we construct a low order Taylor 
polynomial expansion that equals as a point set a shape 
index equivalence class. We begin with an estimate of the 
error that we expect if we expand the signed distance 
function (sdf ) in a 1^^ order gradient Taylor expansion 
plus quadratic remainder term. Let s (x, y) represent the 
sdf over an open neighborhood N{xo/ yo) in the (x, y) • 
plane. 

m 

Then the Taylor esqpansion for s in N(x^, y,) is 



|^Uo.}'o)-(^-^o) + |^(^o.3'o)-(3'-3'o)j + 



2! 



^'■^ iax,ay)-(x-x,y + -ff-iax,ay)-[ix-x,){y-yoy] + 



; a = aix^y) and 0 < Of < 1 



■iax,ay)-{y - yo) 



ay 

We translate s in the z-direction so that s(xo/yo) - 
0, eliminating the constant term. Now rotate the 
orthonormal basis defined by the tangent plane of the sdf 
at (xo, yo) plus the z-axis so that it coincides with the 
orthonormal basis formed the canonical co-ordinate 
system. 

Notes : 
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10 



15 



a- Rotation in the (x,y) coordinate plane will align 
any pair of orthogonal vectors the coordinate 
plane axes. 

b. Rotating the surface so that s(x^,, y^) is a maximum 
causes V sCx^, y©) =0. 

Combining the translation and the rotation, the 
Taylor expansion is 



■ (ax, ay^-x"" + — -5- (ax, ay)'y 



; 0 < <3f < 1, 



dx^' dy' 

The coefficients inside the square brackets are the 
principal curvatures at (cxx,ay). They can be computed 
from the mean and Gaussian curvature at ((Xx,QEy) . 



H (ax, ay) = 



—^(ax,a y) + -—^(ax,ay) 
Bx^ By 



ax" oy 

Using shape analysis, the tme value of every term 
in the Taylor expansion is known. The Taylor disk records 
the first two teims in the expansion, the expansion 

20 point, the disk's radius of convergence, and the maximum 
error incurred on the convergence disk when a sample 
coordinate is approximated by just the constant and 
linear gradient terms. When an intersection curve is 
established, then the definition of the curve is 

25 attached. There is a separate attachment for every 
intersection curve. 
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The a parameter is interesting. Given the center of 
the Taylor disk, the disk radius, and polar angle 
increments then the a parameter describes the 
intersection as well as the principal curvatures along 
5 the intersection curve. (A more descriptive name for this 
parameter is "wireframe".) This parameter determines the 
adequacy of the interpolating fianction in the 3D grid. In 
other words, the adequacy depends on the curvature of the 
parent surfaces along the intersection path. 

10 We prepare the topological analysis of a shape-based 

surface by creating a manifold whose -charts are point 
sets that correspond to shape index intervals. We call 
this a shape index manifold. We explained above how to 
construct the manifold's charts and how to estimate the 

15 differential properties of a chart. These tasks make 
sense with no external context, i.e., a background 
framework. We refer to properties that make sense in this 
self-contained context as "intrinsic". Any property of a 
manifold that is not intrinsic is "extrinsic" . An. exanple 

20 of an extrinsic attribute is the relationship between the 
boundaries of a surface and the boundaries of a shape 
index manifold. 

An intrinsic boundary is one that remains a boundary 
under rigid motion. It is a Cartesian tensor by virtue of 

25 this definition. This bo\andary separates regions that are 
well approximated by patches taken from a surface of 
revolution. An extrinsic bo\indary is a boundary that owes 
its existence to the configuration of backgro\and 
surfaces. Should the background surfaces move then the 

30 shape of the intersection and indeed the very existence 
of the intersection can change. The choice of solver 
depends on the dataset assuitptions . The Shapes ® library 
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does not provide recognition or interaction services 
pertaining to intrinsic boundaries. We agree that a 
characteristic of a soxind approach to point set 
classification is a robust algorithm to compute the 
5 intersection of two extrinsic bovindaries. But we go one 
step further and say that intersection in a region that 
is devoid of extrinsic structure can be still be quite 
complicated, if the intersection involves a non-empty 
subset of a non-dif f erentiable interface that separates 

10 two shapes. 

We seek a data structure that conveys a framework 
overview. We weuit it to contain the complete boundary 
representation for the framework plus a synopsis of the 
shape of every framework surface and curve. We envision 

15 using this data structure to reply to browser level 

framework data base queries when the caller does not want 
to open the framework with the standard geometry services 
toolkit . 

We define a structural synopsis to be a Shapes /GQI 
20 topology graph (also known as a boxindary-representation 
or ""b-rep") plus a shape index manifold description of 
every 2D node in the b-rep.- 

Preferably, a Reeb graph is used to describe a 
configuration's Morse critical points and homotopic 
25 skeleton. See, M. Hilaga, Y. Shinagawa, T. Komura, T. 
Kunii , Topology matching for full automatic similarity 
estimation of 3D , SIGGRAPH 2001, pg. 203-212, and Silvia 
Biasotti, Topological techniques for shape understanding , 
http : //www. eg. tuwien. ac . at/studentwork/CESCG- 
30 2001/SBiasotti/ ■ We recall the definition of a Reeb 

graph. Let M be a path-connected manifold and let f be a 
real-valued on M. Then the Reeb graph associated with (M, 
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f ) is a set of (M, f ) , R) equivalence classes that is 
defined by the relation (x,f{x)) ~ (y,f(y)) if and only 
if f(x) = f{y) and x, y are meiribers of f-1 ( (f (x) ) . In 
practice, the nodes and arcs in a Reeb graph are 
determined from continuous sampling of homotopic 
identification of height field contours. Since the height 
field is a Morse fianction, we obtain infointiation 
•regarding each non- degenerate critical point and the cell 
to which the critical point is attached. The discussion 
of 2D parameterization herein summarizes the relevant 
facts from' Morse Theory. 

Figure 19a shows the Reeb graph of a standard torus. 
Figure 19b schematically illustrates a 2D cell suspension 
that is induced from the axes and planes of symmetry and 
critical point theory. The Reeb graph's nodes in this 
case are critical points on the torus, since the critical 
points are isolated and non-degenerate. The graph 
contains 4 nodes. The bottom and top nodes correspond to 
the height field minimum and maximum. These two points 
lie on an isoparametric curve of constant positive 
Gaussian curvature. The other two nodes correspond to the 
lower and upper saddle points. They too lie on an 
isoparametric curve, but this curve has constant negative 
Gaussian curvature. It is easy to see that the symmetry 
group of the torus leaves invariant the orbit of a saddle 
point as well as the orbit of the min/max. We remark that 
observing this invariant behavior is an easy way to 
discover symmetry transf oinmations . 

A Reeb graph describes homotopic equivalence. Reeb's 
representation has no concept of shape. The Reeb graph 
says nothing about Gaussian curvature or mean curvature 
or shape index. There is no mention of the two opposing 
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curves of zero Gaussian curvature. Nor is there a 
description of the shape of the saddle point curve that 
boxands the interior void, Homotopic eqpjivalence does not 
leave invariant 2D regions of constant shape index, so it 
5 is not possible to reason about symmetry orbits \mder 
Reeb. eciui valence . 

Neither the Hilaga nor Biasotti references augment 
the Reeb graph sTibstrate with curvature information, 
citing the natural instability of curvature estimation in 

10 noise. We agree that curvature estimation in a noisy 

environment is difficult, but we have sibove that it makes 
sense to treat a noisy surface as a smooth substrate that 
is contaminated with noise. We think of the 
representation of a noisy surface as a minimization 

15 problem in Lagrangian mechanics. The smooth approximant 
represents the kinetic energy in the decomposition. We 
define the curvature of the original surface to be that 
of the smoothed component of the original surface. The 
Riemann integral of the discrepancy between the smooth 

20 approximant and the original surface is a measure of the 
potential energy in the original surface. Noise 
introduces imcertainty into the curvature estimate. We 
coittpensate for this uncertainty by working with shape 
index intervals, rather than a point-specific value. We 

25 define a shape synopsis diagram (SSD) to a Shapes /GQI 

boimdary representation where a 2D node is a shape index 
manifold. 

We comment further that Hilaga reports the 
development of a similarity metric for a pair of 
30 triangulated surfaces SI and S2 . Their idea is to create 
a level of detail hierarchy for each surface. The authors 
summarize each level of detail by constructing a Reeb 



48 



0474 GB NP 



graph of the coarsened surface. Hilaga chooses the Reeb 
function to be the integral of geodesic distance measured 
at every vertex of SI and S2 . Hilaga approximates the 
integral using Dijkstra's Algorithm. Hilaga argues that 
this Reeb f\anction is superior to a height field because 
the integral is insensitive to orientation. 

We want to compare Hilaga' s method to the method of 
the present invention. This is not easy. Hilaga works 
with triangulated surfaces rather than implicit surfaces. 
Consequently we must be careful regarding the meaning of 
the term "geodesic path". On a smooth surface we can 
define a geodesic curve to be the straightest possible 
path or the shortest possible path, since the two 
characterizations coincide. On a triangulated surface 
they do not. As is common practice, we select 
"straightest possible" geodesies because an existence 
proof is available for "straightest possible" whereas 
none exists if instead we opt for "shortest possible" 
geodesies. 

We restrict attention to a single surface of 
revolution S, because there exists a theorem that says 
that a parallel on S is a geodesic exactly when the 
tangent vector at any point on the surface is parallel to 
S's axis of symmetry. The symmetry group of S honors this 
constraint, so the geodesic is an orbit under the action 
of the symmetry group. . 

Figure 20 is cross section of the torus shown in 
Figure 19a. Figure 20 shows the toarus opened along a 
planar surface that joins the saddle point circular orbit 
and the min/max circular orbit. 

In Figiire 20, we have overlaid the symmetry group 
orbits associated with critical points. We measure 
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criticality against the standard height field Morse 
function, because geological data is naturally observed 
in depth. We free ourselves of orientation and pose 
limitations by working with the symmetry group orbits of 
5 critical points rather than the isolated critical points 
themselves . 

Figure 21 is a schematic of the shape synopsis 
diagram (SSD) of the torus shown in Figure 19a. The SSD 
reports locations in normal position coordinates and uses 

10 a homogeneous transformation to correctly position this 
data in model space. We also note that an SSD takes 
little storage beyond that already needed to instantiate 
the b-rep. We have placed details of the shape index 
manifold and shape index chart provided below. 

15 Since an SSD augments the standard b-rep we attach a 

spatial frame to the b-rep node to reference the SSD. 
Since the top-level node in the SSD contains a database 
identifier, standard navigation methods can be used to 
locate the b-rep given the SSD. We call attention to the 

20 arrow notation in the SSD. We use an arrow to define the 
correspondence of ID chart boundaries to 2D charts . The 
arrow that is attached to the concave chart signifies 
that the chart does not bound a sub- volume along its 
outside. Similar remarks pertain to the convex chart. 

25 The Reeb graph of a cyclide is identical to the Reeb 

graph of a torus, as expected since the standard Reeb 
graph does not consider shape. The difference between a 
cyclide SSD and a torus SSD is the absence of circular 
orbits connecting critical points. That is, all of the 

30 height field critical points for a cyclide are isolated 
and are fixed by the cyclide 's back to front reflection 
operator . 
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Figure 22 is a cyclide that is shaded according to 
Gaussian and mean curvature. 

Next we construct the SSD for a 2-torus. Figure 23a 
is a bi- torus with its associated REEB diagram. Figure 
5 23b is the visual representation of the bi-torus shape 

synopsis diagram. A shape index manifold is a patchwork, 
with each patch taken from a surface of revolution- So it 
makes sense to display the surfaces of revolution, where 
each surface is decorated with the boiinding curves that 

10 define the patch selection. On the bi-torus 's SSD there 

are two thin circular arcs plus the outer circle identify 
the two toroidal components in the bi-torus. We attached 
the ''slab" label to the connective material between the 
two tori, because we perceive the top of the bi-torus to 

15 be flat, i.e., its Gaussian curvature is zero. The slab's 
internal boundary joins the upper and lower parabolic 
curves on both tori. We have suppressed the ID orbits and 
critical point assignments to keep the diagram readable. 
The bi-torus is a genus 2 sphere, so there is 1 minimum 

20 and 1 maximum. Since the tori are tilted, the height 

field does not see the any critical points along either 
saddle point orbit. Finally, the interior surfaces of the 
slab are concave (hyperbolic) in order to conform to the 
torus ' s exterior convex surface . 

25 The differences between an SSD and a Reeb graph are 

very clear. The use of homotopic identification means 
that the Reeb graph cannot distinguish homotopic figures 
that • have significantly different curvature, e.g., a • 
convex figure from a convex figure. Therefore shape index 

30 analysis is meaningless in the context of a Reeb graph. 

We conclude that although we can map a Reeb graph into an 
SSD, a Reeb graph cannot support a shape 
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index manifold. 

According to a preferred eiribodiment of the 
invention, an efficient hierarchical surface 
5 representation will now be described. We represent an 
implicit surface's signed distance function in a narrow 
band octree encoding of a regularly spaced Cartesian 
grid. We discuss herein how curvature and shape analysis 
of the implicit surface simplify the octree as a data 
10 structure. Smoothing or editing in general is likely to 
change a prior shape analysis, so we will also describe 
herein an adjxinct shape index representation that enables 
fast updating of the octree's information archive. 

15 1. We begin by computing the shape index on the 



implicit surface at the finest resolution in the 
octree. 

2. We partition the shape index data into equivalence 
classes . 



20 



a. We find the finest resolution octants that 
contain the surface's height field critical 



points . 



25 



b. We define the following shape index relation on 
the octant comers that were located in the 
previous step. 



30 



c. We say that two octant comers are related when 
that their respective shape index values are 
contained in the same Cantzler sub division of 
the fiindamental space [-1.0, +1.0] and that 
both values are bounded away from both ends of 



the Cantzler subdivision, 
d. We make a second pass over every equivalence 
class of octree nodes that contains less than a 
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threshold number of sait^les. Let v be such a 
node and suppose that v is on the frontier of 
shape index equivalence classes %! and %2 . 

e. If Xl and x2 are approximated by a single 
elementary shape instance, then we will assign 
the node v to either %! or x2 . 

f . If %1 and X2 are not approximated by a single 
elementary shape instance, then we assign v to 
Xj if the error in the elementary shape 
approximation is less than a threshold 

We introduce the shape index information to the 
coarser levels of the octree. 

a. We say that a low-resolution octant is stable 
if there exists an elementary shape 
representation for all highest resolution level 
nodes that are controlled by the coarse level 
octant . 

b. We record stability status. If a low-resolution 
node is stable, then it will use the elementary 
shape's algebraic evaluators for signed 
distance, mean and Gaussian curvature, etc. 

c. We recognize that a stable low-resolution 
octree node V has the same shape resolution as 
v's finest resolution level components. If we 
want a truly coarser approximation, then we 
develop a separate family of elementary shapes 
for each resolution level in the octree. 

We prune the stable regions of the octree. 
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a- Given an octant Q at level L such that its leaf 
level resolution is part of an ideal elementary 
shape . 

b. We set a status bit in Q that says it coarsely 
5 resolves part of an ideal shape, 

c. Prune the branch of the octree rooted at Q 

d. If Q does not coarsely resolve part of an 
ideal, then we keep the branch rooted at Q. 
Furthermore, we may need to use the leaf- level 

10 interpolator to estimate curvature. 

The situation most favorable to this algorithm is 
when the root surface S* contains a large stable region. 
This can happen if the leaf-level surface S contains a 

15 large region of zero Gaussian curvature, in which case 
the shape index for the region is % or Another 
favorable situation is that both the mean cuarvature and 
Gaussian curvature are constant, e.g., a sphere, so that 
the shape index is again a constant. 

20 Figure 24 is a diagram of the octree with a coarse 

level and leaf level shape index relationship indicated. 
Consistent texture-coding between corresponding regions 
indicates a stable shape index region. 

We consider two examples. The Figure 25 is a diagram 

25 shows part of the French model. A single octree leaf 

defines each plateau. The hemispherical depression has 
constant mean curvature and constant Gaussian curvature, 
so in both regions only one octree leaf is needed. 

Figure 26 is an image of the topography of Crater 

30 Lake, Oregon. Ficrure 26 shows two large plateaus in the 
Southwest part of the lake basin. There are large flat 
regions of the lakebed, so this algorithm will represent 
these regions economically. The small elevation bumps on 
the lakebed will be enclosed in extimded cylinders. 
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We anticipate that the algorithm will perform well 
on the large trimmed conical plateau in the rear. 

According to a preferred embodiment of the 
5 invention, an implicit surface shape identification 

technique will be described. In particular we show how 
to uniformly approximate an implicit surface by a 
patchwork of smooth shapes, where each shape is a section 
of a surface of revolution. The implicit surface is not 

10 required to be smooth. We say that the volume between the 
given implicit surface and the patchwork of smooth shapes 
measures the misfit of the approximation. We improve the 
imiform approximation by reducing the misfit. We reduce 
the misfit by approximating the voliome as a Riemann sum 

15 of generalized prisms. 

Shape identification provides a much higher density 
of information. Figure 27 shows a salt weld in the Gulf 
of Mexico. See, M. Hodgkins, M. O'Brien, Salt sill 
deformation and its implications for subsalt exploration, 

20 The Leading Edge, August 1994, pg. 849-851. We enumerate 
the shapes that collectively represent this complicated 
geobody . 

The shape legend for this image is as follows. 

1 Paraboloid cap 6 

2 . Cylinder 7 

3 Toroidal sector 8 

4 Barrel 9 

5 Conical sector 10 

25 

Now we specify the shape identification algorithm. 
Notation 

30 S A noisy implicit surface. 
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S* A smoothed version of S. 

G A 3D regularly spaced grid 

C A grid cell in G with comers {cO ... c^} . 

a The signed distance fimction defined on C such 

a'^iO) = s n c ^ 0. 
X Tri- linear interpolation associated with grid 

cell C that approximates the signed distance 

function a on C. 
a* The symmetry group of a. 

Shape identification algorithm 



1. We are given a 3D rectangular or rectilinear grid G 
that contains an implicit surface S- 
15 2. Construct a smooth surface S* from eliminating 

noise, but not low freguency intrinsic shape. 
3 • Sample S* such that there exists a sample at each 
entry and exit point on each grid cell face in G. 
We know how to do this, since we know the nature of 
20 a smooth shape that tri-linear interpolation can 

represent . 

4. For each S* sample point estimate the mean and 
Gaussian curvature. 

5. Coitrpute the shape index at every S* sample point. 
25 6- Mark those samples that are close to the center of 

a Cantzler shape index interval. 
7 - Mark those grid cell sairrples that are height field 
critical points, i.e., minima, maxima, and saddle 
points. (We know that these occur at grid cell 
30 comers only.) 

We identify shapes on S*. 

8. For each S* critical point select grid cells 
35 attached to the critical point that have consistent 

mean curvature, Gaussian curvature, and shape index 
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estimators at all sample points in the grid cell 
These grid cell sets define the initial equivalence 
classes . 

9. We want to define a shape from the region of S* 

5 inside the grid cell equivalence class. We check 

that each admitted grid cell's tri-linear 
interpolant has the same symmetry group. If there 
is no common symmetry group, then we siibdivide the 
equivalence class into subclasses that do have a 
10 shared symmetry group. 

10. Using the shape index data cuid the symmetry group, 
identify a surface of revolution that best fits the 
curvature and vinit normal vector field data. 

15 Let C be a grid cell that is iinassigned to a grid 

cell equivalence class, but is adjacent to a grid cell 
that has been assigned to a grid cell equivalence class. 
We want a condition to test for admission to the 
equivalence class. 

20 - . • 

11. There are two categories of grid cells m an 

equivalence class. A member grid cell is either 
"shapely" or is a "misfit" . 

12. In order for the candidate grid cell to be admitted 
25 as "shapely", it iraast satisfy the following 

conditions . 

a. Its tri-linear interpolation function has the 
same symmetry group as all other shapely grid 
cell members. 

30 b. If a shapely grid cell in the equivalence class 

shares a grid cell face with the candidate grid 
cell, then they share a smooth bovindary across 
the grid cell face, 
c. We re-contpute the shape parameter vector and 

35 re-test all grid cells in the equivalence 

class. No "misfit" grid cell is allowed to 
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exceed the misfit tolerance and no shapely grid 
cell can become a misfit. 

13. Test if the candidate grid cell is shapely. 

14. If not, then test if the candidate can be 
5 considered a ^^misfit". 

a. A misfit must be less than some threshold and 
no point on the surface sample in the candidate 
grid cell can be further from the ideal shape 
than some other threshold. 
10 15. If the misfit is acceptably small, then we mark the 

grid cell as a '"misfit" suid admit the grid cell. 

We specify a misfit reduction algorithm later in 
this part of the description. 



15 



20 



It may be that some grid cells have not been 
assigned to an equivalence class. These grid cells are 
not associated with a shape that has a critical point in 
the smoothed implicit surface S* . 



16. Form a new equivalence class by finding connected 
sets of grid cells where the individual grid cell 
tri- linear interpolation functions have the same 
symmetry group and whose shape index samples lie in 
25 the same Cantzler shape index interval. 

It may be that some grid cells can be assigned to 
more than one equivalence class. If we have a choice, 
then we assign the candidate to the class that accepts 
30 the candidate as shapely. 



35 It may be that some grid cells have not been 

unassigned to an existing equivalence class. Since S* is 
smooth, we know that noise is not a contributing factor. 
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Instead, we have loncovered a set of grid cells that 
represent a rapid change in shape. When we encounter a 
grid cell of this type, then we deform the shape in the 
grid cell to conform to any adjacent grid cell and assign 
5 it to the equivalence class that forces the least amount 
of deformation. 

Two ecjuivalence class shapes may fail to connect in- 
a region where they should. When this happens, we deform 
the respective shape in those boundary grid cells tintil 
10 they match the tri-linear interpolation f\inction in each 
grid cell. 

We now discuss the interface between two grid cell 
equivalence classes. If a grid cell face separates the 
two equivalence classes, then we call the interface 

15 "sharp", othein«jrise the interface is "diffuse". If the 
interface is "sharp", then there exists a ID boundary 
separating the two shapes. We compute this boxindary by 
equating the two quadratic surface e3«pressions and 
solving. This is not always smooth, so it will be 

20 necessary to blend a narrow region on each side of the 

intersection curve. If the interface is diffuse, then we 
have a 2D region of overlap between the two shapes in 
which the smoothing can be naturally accommodated. We fit 
a shape to the overlap. We close this topic by checking 

25 this analysis against two adjacent grid cells that have, 
different conic section fibre bundle shape 
approximations. We conclude that the analysis is valid, 
since tri-linear interpolation is continuously 
dif ferentiable on the grid cell face that separates the 

30 two adjacent grid cells. 

Now we overlay the patchwork shape on top of the 
original implicit siirface S and look for chances to 
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reduce the misfit. We will do this by approximating the 
volxime as a 3D Riemaim integral. Our misfit reduction 
procedure is as follows . 

5 Misfit reduction algorithm 

1, V^e project the noisy implicit surface onto 
the patchwork smooth surface. We identify 
connected regions of the implicit surface 
that are responsible for significant misfit 

10 and try to reduce these regions first. 

2. We anticipate multiple misfit regions 
associated to each smooth shape. The most 
complicated element in a misfit description 
is the representation of the outermost ID 

15 boundary of the misfit. We use the outermost 

ID boundary to define a Riemann s\jm 
approximation to the misfit. Consequently, we 
will ignore (initially, at least) low 
aitrplitude but high frequency nearby features. 

20 3. Define height field (relative to the smooth 

approximant) contours at each critical point 
in the noisy intplicit surface. Project these 
contours onto the smooth approximant . 

4. Smooth the set of projected contours- For 
25 each contour in the parameterization we 

identify the critical points on that contour. 
We triangulate . the annular region bounded 
between consecutive contours . 

5. Now we fit a frustum-like collar to the base 
30 of the region where the misfit attaches to 

the background. 

6. The ciomulative shape resembles a telescope, 
hence the name ^"telescopic deformation''. The 
process accommodates tendril-like shapes, in 
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the sense that surfaces cem recursively 
support branches off other earlier branches. 

Figures 28a-c illustrate an exaii^le of the misfit 
reduction process. Suppose that we want to reduce the 
misfit in approximating a tendril shape by a plane, see 
Figure 28a. We subdivide the tendril using height field 
planes that approximately set to the tendril's critical 
points . 

We have generated a sequence of truncated collars 
that reduce the misfit. We are not done, however, since 
this approximation is not dif f erentiable at any point on 
a collar ring separating two truncated cones. We remedy 
this deficiency by replacing a thin swath about the ring 
with a thin swath about the equator of a sphere. 

Figures 29a-c illustrates . an example of blending a 
non-diff erentiable join of two collars. The idea is to 
map given instance (Figure 29a) to the one situation that 
we understand how to fix (Figure 29b) . The one situation 
that we understand how to fix is an isosceles right 
triangle with an inscribed circle such that the center of 
the circle is the centroid of the triangle. We apply the 
linear transf oirmation shown in Figure 29c to the 
isosceles right triangle in panel shown in Figure 29b to 
the exterior and interior triangles - corresponding to Qi 
and 02, respectively - in Figure 29a. 

The misfit reduction procedure is as follows. 

l.We construct height field contours on the noisy 
surface S. Remember that the height field is taken 
normal to the smooth shape approximation S*. 
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2. Working top down, we project each contour to S* . It 
is possible that the projection of successive 
contours cross. 

3. We approximate each contour projection in two steps. 
5 First we locate all ID critical points - minima, 

maxima/ and inflection points. Second we fit a 

smooth approximation to each section of the 

projected contour; where a section is delimited by 

critical points . 
10 4. We triangulate the interior of each projected 

contour, and then construct tetrahedra using 

consecutive contours. 
5 . We approximate the cap /cup corresponding to the 

innermost ring of the bull's eye with the ring's 
15 closure in the plane. 

We next consider how to describe a branching 
structure of hills that separate a region into multiple 
valleys. Figure 30 shows as a geological example a water 
20 breach as indicated by the white arrow. 

Here is a compact description of this branching 
structure. 

1. Find the ridge edge in the branching structure. The 
25 ridge edge is a connected (therefore degenerate) set 

of critical points . 

2 . We assume that as typical cross section is a conic 
section, e.g., a parabola, an ellipse, etc. 

3. We allow for stretching and contracting of the cross 
30 section. We call a region over which the typical 

cross section stretches or contracts a "transition 
zone'' . We provide linear interpolation between the 
start and stop position cross sections. 
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4. We anticipate a valley that runs perpendicularly to 
the ridge edge. For the case of a branching 
structure, we treat a valley of this type as a 
misfit. The following NASA aerial reconnaissance 
5 image contains active breaches in the ridge 

structure. 

We consider a more challenging example. Figure 31 is 
a NASA Shuttle Mission photograph of the Richat Structure 

10 in Mauritania. NASA believes that most likely explanation 
of the origin of this structure is that it is uplifted 
rock sculpted by erosion. NASA also says that there is no 
widely accepted theory that explains why the Richat 
Structure is nearly circular. 

15 The dominant shape in the Richat Structure is a 

nested sequence of toroidal structures. Because of 
erosion, some of the toroidal structures exist as 
sections of a torus only.. Ecker argues on pg. 4-5 of his 
lecture notes that a torus satisfies mean curvature flow 

20 with respect to its mean curvature. This flow is 

singular, in the sense that in finite time the circular 
cross section shrinks to a point and hence the torus 
shrinks to a circle. This is a second example of the fact 
that mcf evolution need not preserve the dimension of a 

25 shape. (The other example is an uncapped cylinder whose 
radius r(t) is given by the formula 

r(f) = ^r^(0)-2f, provided that te (fi,^-^). In finite time the 

uncapped cylinder collapses to its axis of symmetry. , 
With respect to known techniques, Rouby et al- 
30 describe a parameterization algorithm for triangulated 
surfaces. See, D. Rouby, H. Xiao, J. Suppe, 3D 
Restoration of Completely Folded and Faulted Surfaces 
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Using Multiple Unfolding MechanisinS / AAPG Bulletin 84(6), 
June 2000, pg. 805-829. This method involves projecting 
the triangulation onto the (x,y) co-ordinate plane. It is 
not clear how to compensate for triangles that have a 
5 near singular projection. Also, the projection of 

multiple triangles onto a common plane can lead to a 
complicated overlapping and slivers. These are the sort 
of numerical instabilities that we want to avoid and 
which this algorithm does avoid. 

10 

A description will now be give for 3D conformal grid 
generation and reconstruction of shape from disconnected 
pieces, according to the invention. 

Shape identification improves the performance of our 

15 incremental 3D conformal grid generator. Given a sub- 
volume V with boundary 9v = S, we uniformly approximate S 
as a collection of trimmed ideal elementary shapes {Em} . 
Then we enclose S in a tight sequence of rectangular 
prisms {Pk} such that the "vertical" edges of the prism 

20 are normal to the ideal elementary shape approximants . 

We generate a 3D conformal grid in 3 parts - 

1. 3p,^ n S, which is the part between the prism faces 
25 and the measured surface S. 

2 . The bounding box B„ of each trimmed ideal elementary 

shape E„. 

3. A prism P* c V such that P* n B„ = 0^ 

30 We use a proportional spacing correlation scheme to 

grid part #1. In part #2 we use Em's parametric 
expression for normals to grid Bm. In part #3 we grid P* 
in the obvious way. The incremental nature of the 
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procedure follows from the subdivision of the grid into 
regions associated with the prismatic enclosure of S and 
the \iniform approximation by trimmed ideal elementary 
shapes . 

5 Figure 32 shows a conformal grid induced from a 

proportionally spaced correlation scheme. We use a 
proportional spacing correlation scheme to define the 
grid iso-surf aces . Here is an ideal application of this 
mechanism (This is an idealization of a cross section 

10 through the Gullfaks field.). We see that the sub-blocks 
have approximately the same z-direction thickness. It is 
irrelevant how the horizontal extents correspond- In 
other words, we can generate either a regular spacing 
grid or an irregular spacing grid. We remark that a 

15 compromise such as a '^tartan grid" .{also known as a 

rectilinear grid) is more memory efficient than a regular 
grid but still has a convenient algebraic lookup 
function. 

We observe that a correlation scheme implements a 
20 cheap form of mean cuirvature flow and if generated 

indepenjiently of a parameterization, then the conformal 
grid is a roundabout way to also generate a 
parameterization . 

25 We conclude this part of the description with a 

discussion of a way to convert an existing 3D Cartesian 
grid to a conformal grid. Figure 33 illustrates a non- 
conformal 3D Cartesian grid. As is well known, it is non- 
trivial to trim grid cells that cross volume boundaries. 

30 Here the backgrovmd grid shown in dashed lines, crosses 
the overburden. It is clear in this diagram that the 
background grid cells do not conform to lithological 
boundaries . 
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We remedy these deficiencies in the following way. 

!• We attach the volume of interest to a spatial 
frame . 

5 2 . We conf onnally grid the voliame that lay between 

the overburden and the reservoir. 
3, For best results we want a transition zone from a 
Cartesian grid to the conformal grid, preventing 
reverberation in the grid across lithological 
10 boundaries. 

Forensic reconstruction will now be discusses. 
Suppose that an application identifies a set of 
disconnected 2D patches and would like to construct a 

15 surface from the patch set. We show how mean curvature 
flow enables us to construct such a surface - hence the 
sub- title for this discussion. This techniq[ue will be 
seen to be similar in execution to grid generation, which 
is why we discuss it here. 

20 Suppose that we have a collection of 2D patches in 

R^. We assxime that all of the patches are oriented in a 
consistent manner, which must be the case if all of the 
patches are taken from a common surface. We construct the 
solution using patches or parts of patches that have a 

25 common Gaussian curvature signum. We do not insist that 
the Gaussian curvature signum be constant across the 
patch point set. However, if it is not constant, then we 
use only those patches and parts of patches that have the 
same signum to construct a surface. This procedure can 

30 create holes in the surface. We fill in these holes with 
the remainder of patches that were used only partially. 
Some patches with the opposite signum may not fill in 
holes, but instead form folds in the surface. We 
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construct a second surface from these patch point sets 
and union the two parts together. 

We remark that the following algorithm is applicable 
to the task of sewing together parts of the bounding 
5 surface of a geobody, e.g., a salt body. 

Our method is as follows . 

1. We assume that the velocity of the flow equals l,so 
10 that a time interval is equal to a distance. 

2. For each patch or region of a patch R that has ( + ) 

signtam Gaussian curvature, we define a uniformly 
e3«panding mean curvature flow with initial 
condition equal to R. 
15 3. These patches are taken from a common surface, so we 

saiiple each flow at the same time step increments. 

4 . we stop a particular flow when it intersects another 

flow pattern. 

5. Eventually each flow intersects a different flow. 
20 6. It may be that the confluence of two flows is not 

dif ferentiable. If this happens, then we define a 
flow to smooth the sharp angle. Alternatively, we 
can smooth the sharp angle using the Squeeze 
operator defined in the Misfit Reduction 
25 algorithm. 

If a formation contains a t\annel, then we may wish 
to reverse the dissolution process, i.e., fill in the 
tiinnel. In this situation we will approximate the timnel 

30 as a collection of uncapped cylindrical and toroidal 

sections. Ecker's lecture notes show that mean curvature 
flow applied to an uncapped cylinder will shrink the 3D 
surface onto its ID axis of symmetry in finite time. 
Smocyzk has obtained a comparable result for a torus. His 

35 result shows that a torus collapses vmder mean curvature 
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flow to its inner parabolic circle in finite time. We 
conclude that when applied to a tiinnel, mean curvature 
flow will collapse the tunnel to a curve in finite time. 
•Figure 34 is an image of the Devil's Potholes , South 
5 Africa. We can identify two different kinds of pothole 
formations in this image. The circled pothole on the 
right is well described as a cylinder. We can collapse it 
to its axis of symmetry - which is a line - under mean 
curvature flow. The pothole on the left is more complex. 

10 The top of this region is a curved throat, which is 
followed by a sharp edged cone. Finally there is a 
toroidal section that is turning away from the viewer. We 
know that each section will collapse to a singular edge 
in a finite amount of time. Taking the longest time 

15 interval needed to collapse a section of this structure, 
we guarantee that the entire pothole will collapse to a 
polygonal path. 

We now discuss a related problem. Given an image of 
a layered sequence, suppose that an unconformity. The 

20 unconformity will erase everything in the image that is 
above the \anconf ormity . How good a reconstruction of the 
eroded interface can I achieve, given just the migrated 
image. As an example, we return to the diagram of a cross 
section of a shallow sequence in Turkmanistan shown in 

25 Figure 3. We are interested in the guess of what the 
missing section of the layer looks like. 

Our algorithm assumes that the formation satisfies 
the following assumptions. . 

30 1. There exists a reference layer such that the visible 
part of the sequence is well approximated by a 
single reference surface adaptively sampled 
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correlation scheme. In the above schematic, the 
vertically striped layer R is the reference layer, 



2. We reconstruct the eroded region of a layer by 
5 evolving the reference interface under mean 

curvature flow until the front joins the 
\anconfonnity. In the above schematic, the 
reconstructed interfaces are shown as dotted 
outlines above the \inconf ormity. 

10 

While it may not be obvious from Figure 3, conformal 
grid generation is a very short time solution to mean 
curvature flow. In other words, if erosion had not been 
iit^>osed then we would recognize the conformal grid 

15 generation in the framework close to the reference 

surface and flattening further away. As an illustration, 
we consider Figure 35, which is a NASA LANDSAT image of 
the Yukon River delta (NASA Geomorphology plate D-9) . We 
focus on the boundary cracks that subdivide the roxinded 

20 joint at the left top of the image. These boundaries 

delimit curvature flow regions that collectively form the 
bend in the deltaic mass. It has been proposed to use a 
diffusive process. See, J. Davis, S. Marschner, M. Garr, 
M. Levoy, Filling holes in complex surfaces using 

25 volumetric diffusion . First International Symposium on 3D 
Data Processing, Visualization, and Transmission, Padua, 
Italy, June 19-21, 2002. This reference does not make 
use of curvature information in the algorithm used. 
Furthermore, the reference does not deal with closing 3D 

30 voids such as tunnel closure. 
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We have discussed how to construct a patchwork 
covering of an iinplicit surface from, first order Taylor 
polynomials. The error term is cjuadratic with 
coefficients that are the principal curvatures defined at 
5 a point that is somewhere within the circle of 

convergence about the expansion point. The radius of the 
circle is chosen so that the error term is smaller than a 
user-defined tolerance. The error estimate is not sharp, 
because it is not obvious where to evaluate the error 

10 term's principal curvature values, i.e., what is the 

definition of a. Since this error estimate is not sharp, 
we supplement the Taylor representation with a fast 
technique to evaluate a. 

From curvature analysis we derive a representation 

15 of a surface as a connected sum of trimmed elementary 

shape. Some examiples of elementary shapes are (sections 
of) an ellipsoid, hyperboloid, cyclide, and a prism. An 
elementary shape has algebraic expressions for mean 
curvature and Gaussian curvature. These expressions can 

20 be sxabstituted into the quadratic formula to generate an 
algebraic expression for principal curvature. 
Unfortunately the resulting expression is often so 
complex algebraically that a symbolic algebra package 
such as Mathematica® must be used for any manipulation 

25 other than simple evaluation. This is awkward in an 

Engineering environment, so a low order Taylor series 
approximation is preferable. 

Here is an intersection algorithm for implicit 
surfaces. This algorithm assumes that every surface's 

30 signed distance f\inction has been defined on a shared 
grid. 
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Notation 
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15 

Intersection algorithm 

1. Find all grid cells {CJ that surface intersects - 
A grid cell is involved if and only if the signed 

20 distance function is not constant at all grid cell 

corners . 

2. In this algorithm we want to work with surface 
patches that project invertibly onto a face of grid 
cell. We siabdivide the part of a surface contained 

25 in a grid cell into the desired patchwork by 

checking the surface normal field restricted to the 
grid cell. We place a point on the surface in a 
patch exactly when the xanit normal at the point is 
within tolerance of being parallel to all other 

30 surface points that are already assigned to the 

patch . 

For each patch find all S, patches such that their 
corresponding sdf patterns do not rule out intersection. 
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For clarity we assiame that each surface contains exactly 
one such patch. 

ZjOOP$0: 



3. Let Co be the midpoint of the line segment 

connecting two faces of C such that the signuin of 
the sdf for on one of the connected faces is 
different from that on the other. 
10 4. Let j = 0. 

IaOop$l: 

5. Project c^ onto and 5^, respectively. 
15 6. Label the projection s^^ and s^^, respectively. 

7. Either s^^ = Sj^ or not. 

8. If s^j = Sjj, then go to Lal>el$2. 

9 . If s^j ^ Sjj^ then construct the line segment X^^^ in C,, 
connecting s^j to s^j . 

20 10- Let Cj^i be the midpoint of X^^^ 

11. It might be the case that c^ « Cj^^. We expect to see 
this when a very thin hyperbolic neck is contained 
in the same cell as C^. 

12 . If Cj ~ c^^^, then we want to jump away from this 

25 liopeless local minimum. We do this by comparing the 

signum of the sdf for Sj^^ and the signiam of the sdf 
for Sjj to signum of sdf obtained on the grid cell 
faces joined by the line segment X^. Moving in the 
direction that differs from the current values, jump 

30 to the midpoint along A^^,. 

13. Increment j and return to Loop$l. 

I.al>el$2: 
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We use a Newton algorithm to find another point of 
intersection^ 

Loop$3: 

14. Let (x^rYii) a point of intersection for and S, 
in C and suppose that (x^*,y^*) is another point of 
intersection in C. Expanding each surface in a l" 
order Taylor series about (x^/Yia) / we have 

52U12*. 3'i2*) = '^2(^12. 3'i2) + V52(^2. ^'12) • (^12 - ^12*' - ^12*) 

15. Siibtracting, we have 2 linear equations in two 
\anknowns. Solve for the non-trivial solution 

0 = V - 5i )(j:i2, jiz) • (^2 - ^12*' yi2 - y^*) 



16. Return to Loop$3 ^mtil (x,,*, y,,*) is acceptably 

close to being another point of intersection in C. 
20 17. We compute the intersection cuirve between (.x.^, y^) 

and (x„*, y^j*) by repeating the midpoint svibdi vision 
algorithm in I<oop$2. 

18. Return to Loop§3 and repeat the l" order Taylor 
series algorithm in Loop$2 on both ends {x^,yi2) and 

25 (Xij*,yij*) xintil the curve crosses the faces of C^. 

19. Return to Loop$3and continue processing grid cells 
in {C^} vmtil exiiaustion. 



Visualization of an intersection curve is as follows. 

1. Project the intersection curve onto the internal 
cylindrical parameter space. 

2. Facet the cylinder using a constrained Delaiinay 
algorithm. 
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3. Invert the mapping, projecting the triangulation of 
the region enclosing the ^'f lattened'' intersection 
curve to the parent implicit surface. 



placing sample points on a surface and generating a 
narrow band octree from them. Now we discuss how spatial 
frames can be used to implement an efficient update 
mechanism. For details of the reappearance of SIGMA' s 
10 hybrid grid-mesh concept, refer to U.S. Patent 

Application Serial No. 09/163,075, incorporated herein by 
reference . 

Each geological feature has its own octree which 
overlays part of the background octree. An overlay can 

15 cover a strictly smaller region of the background, and 

this small region can be sampled differently (adaptively) 
from the backgrovmd. The sampling can be adapted to the 
coinplexity of the framework region that the overlay 
encloses. The overlay changes as the framework region 

20 changes. A separate spatial frame is assigned to control 
each octree overlay- Overlays are loaded on demand. 
Within a spatial frame the sampling can be adjusted in an 
adaptive manner. A frame boundary sample node can be part 
of multiple octree overlay. An entity that is outside a 

25 particular spatial frame learns about the frame's 

enclosed shape by interrogating the enclosing spatial 
frame. Given an octree identifier and a location inside 
or on the botindary of .the octree's enclosing frame, the 
local cpt returns the cpt and sdf for the remote octree 

30 overlay. Adjacent frames might support a different set of 
ajad number of framework surfaces. 

Given a region of interest, we associate to each 
spatial frame the point set that the frame contains. We 
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We have described a curvature-adaptive method for 
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represent the set of spatial frames associated with a 
region of interest as a directed acyclic graph. Two 
parent frames that partially overlap will be represented 
in the graph as nodes that point to a common descendant. 
Thus spatial frames constitute a topology graph. 

As an example, a vsp was acquired over part of the 
Zechstein Salt formation. Figures 36a and 36b show, for 
reference, the background Zechstein Salt and the region 
in the Zechstein where the vsp was acquired. Figure 37 
show the frame graph that ties the vsp region of interest 
to the Zechstein Salt backgroiind. We use the inner frame 
to outer frame relationship to express «A is a boundary 
of B" and "A is logically part of B" . In GQI terminology 
a spatial frame as defined herein is equivalent to a 
15 gqi_Frame_t. We have extended this fxanctionality two 
ways. First, we maintain logical containment, i.e., 
feature relationships. Second, we allow mixing space and 
time frames . 

Spatial frames are useful in siabdividing a faulted 
region to match regions svibject to different physics. For 
exanple. Figures 38a, 38b, 39a and 39b illustrate the 
separation of faulted sediments from unfaulted sediments. 
Figure 38a and 38b show the case of a sequence of. 
sediments such that some but not all of the sediments 
have been faulted, and a simple normal fault is involved. 
Figures 39a and 39b show the case of a fault network 
where some of the faults in the network emanate from 
other faults. Figure 39b shows that the \xnique final 
configuration of spatial frames that partition the 
voliames shown in Figure 39a. Using mathematical 
induction on the nuiriber of faults in a compact region of 
interest, we conclude that it is immaterial in what order 
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spatial frames are assigned to isolate distinct 
sedimentary regions. In other words, the final 
subdivision is always the same no matter what the order 
of isolation is. 
5 We think of mean cuirvature flow as an evolutionary 

process. For siittplicity, we assume that any evolutionary 
process is a structure group of dif f eomorphisms of some 
4D fibre bundle. All information regarding the process at 
a moment in time is encoded in a time frame. A time frame 

10 contains a reference to a region of interest that in turn 
is represented as a 3D fibre bundle. We do not demand 
that an evolutionary process provide a physically 
plausible explanation of a formation state, rather it is 
enough if the visual expression appears plausible. 

15 Figure 40 illustrates a time-lapse seismic evolution 

(steam injection front tracking, with the left panel 
^'bef ore'' and the right panel ^^af ter'' . ) . We discuss now 
how we use spatial frames to define a topology graph. Let 
S be an implicit surface and let {Ek} be a family of 

20 elementary shapes that approximate S. Each E^ reduces its 
misfit with S by fitting a disjoint region of S to a 
telescopic extrusion of E^. - 

We use a different spatial frame to contain each 
telescopic extrusion and attach the collection of these 

25 frames. We assign each elementary shape to its own 

spatial frame and refer to the misfit reduction through 
attachment of the spatial frame list to the parent 
elementary shape frame. Figure 41 shows the reference 
structure for this set of relationships. 

30 

We summarize the contents of the new class instances 
that we have defined in this paper. 
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